//============================================================================
// Name        : csiena
// Author      : Josh Lospinoso
// Copyright   : MIT License
// Description : Austere, limited implementation of the R package RSiena in
//				 C++. Used for fast maximum likelihood and method of moments
//				 estimation of stochastic actor oriented models (SAOMs)
//				 without any frills or user experience niceties.
//============================================================================

#include "globals.h"
#include "chain.h"
#include "mersenneTwister.h"
#include "mcmcmle.h"

using namespace std;
using namespace csiena;

void simCompleteDataEstimation(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(50, atoi(argv[1]));

	e->loadStart("data/s501.csv");

	cout <<"b_od,s_od,i_od,b_re,s_re,i_re,b_tt,s_tt,i_tt,b_3c,s_3c,i_3c" << endl;
	for (int iter = 0; iter < atoi(argv[2]); iter++) {
		e->setBeta(0, -1.0);
		e->setBeta(1, 1.0);
		e->setBeta(2, 0.25);
		e->setBeta(3, -0.30);
		e->generateSimulatedChain(250);
		for (int i=0; i<500; i++) {
			e->estimateCompleteData();
			e->partialSingleCoordinateNewtonRaphsonStep(.01);
			//e->partialNewtonRaphsonStep(.01);
		}
		printf("%f,%f,%f", e->getBeta(0), e->getBetaScore(0), e->getBetaInformation(0,0));
		for (int j=1; j<BETA_LENGTH; j++) {
			printf(",%f,%f,%f", e->getBeta(j), e->getBetaScore(j), e->getBetaInformation(j,j));
		}
		printf("\n");
	}
}


void simEstimation(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(50, atoi(argv[1]));

	e->setBeta(0, -2.47);
	e->setBeta(1, 1.40);
	e->setBeta(2, 0.95);
	e->setBeta(3, 0.75);
	e->loadStart("data/s501.csv");
	e->generateSimulatedChain(150);
	e->completeData->print();
	e->executePhaseOne(true,false,true,2500,0,1000, .1, true);
	e->executePhaseTwoSubphase(25,250,.06,true,0,100, true);
	e->executePhaseTwoSubphase(25,350,.03,true,0,100, true);
	e->executePhaseTwoSubphase(25,450,.015,true,0,100, true);
	e->executePhaseThree(0,5000,true);
}

void s50estimation(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(50, atoi(argv[1]));
	e->setBeta(0,  -2.30); // Density
	e->setBeta(1,  2.0); // Reciprocity
	e->setBeta(2,  .50); // Transitive triplets
	e->setBeta(3,  0.0); // In degree popularity
	e->setGamma(0,.01); // False negative
	e->setGamma(1,.025); // False positive

	e->loadObservation("data/s501.csv");
	e->loadObservation("data/s502.csv");
	e->loadObservation("data/s503.csv");
	cout << "Loaded s50 dataset." << endl;
	e->completeData->evaluateChain();
	e->completeData->print();

	e->executePhaseOne(true,true,true,500,0,1500, 1, true);
	/*
	e->executePhaseTwoSubphase(300,550,.2,true,0, 500, true);
	e->executePhaseTwoSubphase(300,500,.1,true,0, 500, true);
	e->executePhaseTwoSubphase(300,450,.05,true,0, 500, true);
	e->executePhaseTwoSubphase(300,450,.025,true,0, 500, true);
	e->executePhaseThree(0,3000,true);
	*/
}

int simulateScoresFromS50(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(50, atoi(argv[1]));

	e->setBeta(0, -1.10);
	e->setBeta(1, 1.1);
	e->setBeta(2, 0.3);

	e->loadStart("data/s501.csv");
	e->generateSimulatedChain(200);

	// For simulating from s501 and gathering scores, information.
	// Number of simulations:
	int iterations = atoi(argv[2]);
	// Number of score estimates to be taken per iteration:
	int samples = 1;
	int mcmcSamples = 500;
	int mcmcSpacing = 0;

	e->loadStart("./data/s501.csv");

	// Simulate an observation s502* from the s501 dataset, and collect scores/info from conditional draws
	// from intervening chains via MCMC.
	cout << "OScore,RScore,TScore,CScore,OCInfo,RCInfo,TCInfo,CCInfo,OOInfo,ROInfo,TOInfo,COInfo,AR,Iteration" <<endl;
	for (int i=0; i<iterations; i++) {
		e->generateSimulatedChain(150);
		//e->burnChain(2500);
		for (int n=0; n<samples; n++) {
			e->estimateMCMC(mcmcSpacing, mcmcSamples, true);
			for (int k=0; k<BETA_LENGTH; k++) {
				cout << e->getBetaScore(k) << ",";
			}
			for (int k=0; k<BETA_LENGTH; k++) {
				cout << e->getBetaInformation(k,k) << ",";
			}
			for (int k=0; k<BETA_LENGTH; k++) {
				cout << e->getObservedInformation(k,k) << ",";
			}
			cout << e->getProposalAcceptanceRatio() << "," << i << endl;
		}
	}
	return 0;
}

int profileLogLikelihood(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(50, atoi(argv[1]));

	e->setBeta(0, -3.10);
	e->setBeta(1, 1.7);
	e->setBeta(2, 1.0);

	e->loadStart("data/s501.csv");
	e->loadObservation("data/s502.csv");
	e->burnChain(10000);
	int samples = 100;
	double pll[samples];
	double b[samples];
	e->estimateProfileLogLikelihood(2,-3,5,samples,b,pll, 100, 10);
	for(int i=0; i<samples; i++) {
		cout << b[i] << "," << pll[i] << endl;
	}
	return 0;
}

int checkChainsCompleteData(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(50, atoi(argv[1]));
	e->setBeta(0, -2.10);
	e->setBeta(1, 2.0);
	e->setBeta(2, 0.05);
	e->setGamma(0, .05);
	e->generateSimulatedChain(200, 150, 50);
	e->completeData->print();
	e->estimateCompleteData();
	e->print(true);

	return 0;
}

int checkChainsMCMC(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(50, atoi(argv[1]));
	e->setBeta(0, -1.0);
	e->setBeta(1, 1.0);
	e->setBeta(2, 0.15);
//	e->setGamma(0,.05);
	e->setAlpha(6.0);
	e->generateSimulatedChain(200,0,3);
	e->executePhaseThree(0,100,true);
	return 0;
}

int checkCopyChains(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(10, atoi(argv[1]));
	e->setBeta(0, -1.10);
	e->setBeta(1, 1.0);
	e->setBeta(2, 0.5);
	e->setBeta(3, -0.3);
	e->generateSimulatedChain(100);
	e->completeData->print();
	Chain* ichain = new Chain(*(e->completeData));
	e->burnChain(75);
	cout << "Printing out the copy:" << endl;
	ichain->print();
	cout << "Reevaluating the chain and printing:" << endl;
	ichain->evaluateChain();
	ichain->print();
	cout << "Finally compare this with the chain we burned" << endl;
	e->completeData->print();
	delete(ichain);
	return 0;
}

int checkMultithread(int argc, const char* argv[]) {
	mcmcmle *e = new mcmcmle(50, atoi(argv[1]));
	e->setBeta(0, -1.10);
	e->setBeta(1, 1.0);
	e->setBeta(2, 0.5);
	//e->setBeta(3, -0.3);
	e->generateSimulatedChain(250);
	e->estimateMCMC(0,1000, true);
	e->print(true);
	return 0;
}

int main(int argc, const char* argv[]) {
	s50estimation(argc, argv);
	return 0;
}

